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EFFECT  OF  THE  TEMPERATURE  OF  THE  MODERATOR  ON  THE  VELOCITY  DISTRIBUTION 
OF  NEUTRONS  WITH  NUMERICAL  CALCULATIONS  FOR  H  AS  MODERATOR 


By  E.  P.  Wigner  and  J.  E.  Wilkins,  Jr. 


ABSTRACT 

In  this  paper  we  set  up  an  integral  equation  governing  the  energy  distribution  of  neutrons  that 
are  being  slowed  down  uniformly  throughout  the  entire  space  by  a  uniformly  distributed  moderator 
whose  atoms  are  in  motion  with  a  Maxwellian  distribution  of  velocities.  The  effects  of  chemical  bind¬ 
ing  and  crystal  reflection  are  ignored.  When  the  moderator  is  hydrogen,  the  integral  equation  is  re¬ 
duced  to  a  differential  equation  and  solved  by  numerical  methods.  In  this  manner  we  obtain  a  refine¬ 
ment  of  the  dv/v3  law. 

*  *  *  *  * 


1.  There  are  two  particularly  simple  problems  in  connection  with  the  energy  distribution  of 
neutrons  that  are  present  in  a  medium  of  finite  temperature.  In  the  first  problem  the  slowing  down 
is  uniform  throughout  the  entire  space  that  is  itself  uniformly  filled  with  the  slowing  down  material. 

In  this  case  the  neutron  distribution  is  evidently  the  same  over  all  space.  In  the  second  problem  the 
neutrons  enter  a  half  space  from  one  side  with  uniform  intensity  and  diffuse  into  it.  The  question  in 
this  case  is  the  density  distribution  of  neutrons  at  large  distances  from  the  boundary  plane  of  the 
half  space  and  the  exponential  relaxation  length  of  the  neutron  density.  We  shall  be  interested  only 
in  the  first  problem. 

There  are  three  phenomena  which  make  the  calculation  of  the  energy  distribution  of  the  neutrons 
somewhat  complicated.  These  are  the  finite  velocity  of  the  particles  with  which  they  collide,  the 
effect  of  the  chemical  binding  on  the  scattering  cross  section,  and  finally  the  effects  of  crystal  re¬ 
flection.  We  shall  disregard  the  last  two  effects,  and  the  calculation  to  be  given  will,  therefore,  be 
valid  only  in  a  monoatomic  gas.  Seitz  and  Goldberger  are  engaged  in  the  study  of  the  crystal  effect, 
and  Teller  has  made  considerations  on  the  effect  of  the  chemical  binding. 

To  derive  an  equation  satisfied  by  N(v),  where  N(v)dv  is  the  number  of  neutrons  per  cm3  whose 
velocity  lies  between  v  and  v  +  dv,  we  proceed  as  follows.  Let  us  denote  the  probability  that  a  neutron 
with  velocity  v^  will  acquire  in  unit  length  of  time  a  velocity  between  v  and  v  +  dv  by  collision  with  an 
atom,  by  NasP(v,v^)  dv.  Here  N  is  the  number  of  atoms  per  cm3  and  os  is  the  scattering  cross 
section  of  the  moderator.  Hence  the  number  of  neutrons  acquiring  a  velocity  between  v  and  v  +  dv  in 
unit  length  of  time  is 

Ndv  jf  ffgFCv.Vj)  N(vx)  dvr 

Let  N  orsV(v)  be  the  probability  that  a  neutron  with  velocity  v  will  be  scattered  in  unit  time,  and  let 
N  asy  be  the  probability  that  a  neutron  with  velocity  v  will  be  absorbed  in  unit  time.  If  we  adopt  the 
1/v  law,  this  second  probability  will  be  independent  of  the  velocity  of  the  neutron.  We  further  assume 
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that  the  scattering  cross  section  is  independent  of  the  velocity  of  the  neutron.  This  is  a  valid  assump¬ 
tion  for  a  moderator  gas,  since  our  results  will  be  applied  only  in  a  reasonably  close  neighborhood  of 
thermal  energies.  It  is  now  clear  that  N(v)  satisfies  the  following  integral  equation: 


f  P(v>vx)  N(vj)  dvx  =  [ V(v)  +y]  N(v). 


Before  calculating  P  and  V  we  shall  derive  a  property  of  the  above  equation  that  simplifies  cal¬ 
culations  greatly  and  also  can  serve  as  a  check.  Evidently  equation  (1)  is  satisfied  by  the  Maxwell 
distribution  if  y  =  0.  This  together  with  the  principle  of  detailed  balance  permits  us  to  give  equation 
(1)  a  symmetric  form.  If  we  denote  the  Maxwell  distribution  of  velocity  by  M(v) 

4$  3  n  2  2  n  i 

M(v)  dv  =  —  v2  e~&  v  dv,  jT  =  (2) 

yr!T 


the  principle  of  detailed  balance  reads 


P(v,vx)  M(vx)  =  P(vx,v)  M(v). 


Herein  the  mass  of  the  neutron  is  used  as  the  unit  of  mass.  It  follows  that 


-  1  P(v,v,)  V M(v.)  =  1  P(v,,v)  VM(v) 

v  M(v)  V  M(vj) 

1  (3a) 

=  =  S^.v) 

is  a  symmetric  kernel.  Thus  if  we  introduce  the  new  quantity  v(w)  =  N(v)/VM(v)  into  the  equation  (1), 
we  get  the  equation 


/s(v)Vl)  r(vx)  dvx  = 


[  V(v)  +  y]  v(v) 


which  has  a  symmetric  kernel. 

2.  The  next  task  is  to  calculate  V  and  the  kernel  P.  Both  quantities  depend  of  course  on  the 
temperature  of  the  gas  and  also  on  the  mass  of  its  atoms.  The  calculation,  although  quite  laborious, 
is  entirely  straightforward  and  gives  for  v^  <  v 


P(v,v1) 


_  26  v  2  _  „2 


v 


e^  (v-^-v2)  [l(pevx  -P(v)  +  l(pev1  +  /3fv)] 


Y~n  v1 


[l(/36>v -/3fv1)  -I(j3f vx  +  00v)]  . 


For  v^  >  v  it  gives 


P(v,vj) 


yr~7T  Vj^ 


-  [l(/30v  -/3f  Vj)  +  l(B9v  +  /SfVj)] 

iT*  vx 

■  e^2(vi2-v2)  [i^evj-^fyJ-i^fv  +  jSev^] 


(5a) 


AECD  -  2275 


[3 


Herein 


,  p- 


2ym  2  \TET 

where  m  is  the  mass  of  the  atoms  of  the  moderator  in  units  of  the  neutron  mass.  I  is  the  odd  function 


(X)=  / 


e~x  dx. 


For  V  we  obtain* 

2  2 

V(v)  =  —  (2v  +  1/02  mv)  1(0  l/inv)  +  e'^  mv  /0  (6) 

yHT 

Equations  (5)  and  (6)  determine  the  quantities  which  occur  in  the  integral  equation  (1).  In  the  deriva¬ 
tion  it  was  assumed  that  the  scattering  cross  section  as  is  independent  of  the  velocity  of  the  neutrons. 
This  is  a  valid  assumption  for  a  moderator  gas  since  (1)  will  be  applied  only  in  a  reasonably  close 
neighborhood  of  thermal  energies.  It  was  further  assumed  that  the  scattering  is  spherically  sym¬ 
metric  in  the  center  of  mass  coordinate  system.  The  physical  significance  of  y  is  y  -  o^v/  org  where 
a_  is  the  absorption  cross  section  for  the  relative  velocity  v  between  neutron  and  atom — y  is  in- 
dependent  of  the  latter. 

The  equations  (5)  and  (6)  can  be  simplified  considerably  by  measuring  the  velocity  in  terms  of 
the  thermal  velocity.  Let  us  denote  therefore 

0v  =  x,  0v,  =  x,  (7) 


r  =  0  y 


orsJ/2kT 


Then  r  is  the  probability  of  absorption  of  the  neutron  while  it  goes  one  mean  free  path  with  the 
velocity  corresponding  to  the  energy  kT. 


We  then  have 


P(x,x1)  =  - —  e^Xl  X^  [l(0x,  -  £x)  +  l(0x,  +  £  x)] 

-  1  VT  xi 


282  x 


+ -  -  [l(0x  -  (  xi) -l(( x*  +  0x)]  for  x*  <  x 

w  xi  1 


2  e2  x 


P(x,Xj)  = - —  [l(0x  -  £  xt)  +  l(0x  +  £ xj)] 

yir  i 
yir  i 


*  V  is  also  given  by  G.  Jager  in  Winkelmann’s  Handbuch  der  Physik  (Barth,  1906),  vol.  3,  p.  698. 
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for  >  x.  We  further  hav* 


and  the  integral  equation  becomes 


—  (2x  +  l/mx)I(m2x)  +e~mX  /Y  m u 


J  S(x,x1)  Kxp  cb^  =  [  V(x)  +  r]  v  (x). 


Although  the  derivation  of  the  expressions  for  P  and  V  is  quite  straightforward,  it  may  be  worth¬ 
while  to  give  a  few  of  the  intermediate  steps  to  facilitate  checking  the  equations. 

V(v)  is  the  probability  of  a  collision  of  the  neutron  with  an  atom  of  unit  cross  section  and  the 
velocity  distribution 

Mm(v2)  dv2  =  4m  -Z  v22  e_/3  ^  dv2-  (10) 

YT 

The  probability  for. a  collision  with  an  atom  of  velocity  v2  is,  if  the  neutron  velocity  is  v,  and  the 
cosine  of  the  angle  between  the  directions  of  motion  is  p,  simply 


vr  =  V  y2  +  v22  “  2w2/i 


which  is  the  relative  velocity  of  the  tw7o  particles.  The  number  of  atoms  within  unit  velocity  range  at 
v2  is  given  by  (10),  the  probability  that  }i  lies  between  \i  and  p  +  dp  isidp  so  that 

1 


V(v)  =  J  |d \i  J  dv2  Mm(v2)  j/v2  +  v22  -2vv2p 


The  integration  over  (i  is  elementary  and  gives 


/  [(v  +  v2>3-|v-v2|3]dv2- 


When  integrating  this  expression  one  must  proceed  with  the  two  cases  v2  <  v  and  v2  >  v  separately. 
This  kind  of  disjunction  is  characteristic  also  for  the  calculation  of  P.  With  this  proviso,  the  integra¬ 
tion  of  (10c)  can  be  carried  out  easily,  and  it  gives  (6). 

The  calculation  of  P  is  more  cumbersome.  Collision  of  a  neutron  with  the  velocity  v^  with  an 
atom  of  velocity  v2  gives  a  velocity  between  v  and  v  +  dv  for  the  neutron  with  the  probability 


P(v;  vx,  v2,  p)  =  0 


v  2  —  v  -2 
vmax  vmm 


if  v  <  vmin 
d  vmin  <  v  <  vmax 


if  v  >  v 


These  equations  express  the  familiar  fact  that,  after  the  collision,  the  probability  is  uniform  in  energy 
scale  between  the  minimum  and  maximum  energies.  Herein,  p  is  again  the  cosine  between  the  direc¬ 
tions  of  motion  of  neutron  and  atom  before  the  collision, 
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vr  -  mvr/(m  +  1), 


vmax  =  vc  +  mvr/(m  + 


c  m  +  1 


Vj2  +  m2v22  +  2mv1V2ju 


is  the  velocity  of  the  center  of  mass  of  neutron  plus  atom,  mvr/(m  +  1)  is  the  velocity  of  the  neutron 
in  the  center  of  mass  coordinate  system.  One  derives  (11)  most  easily  by  a  geometrical  argument. 


Pfv.v^  / idp  J  dv2  Mm(v2)vr  P(v;  v1;  v2,  m). 


=  (  1  4  ”)2V  //  dv2  Mm<v2)  (yl2  +  “S2  +  2mW>' 


where  the  integration  is  to  be  extended  over  the  region  limited  by  -  1  <  \x  <  1  and  that  vmin<  v  <  vmaY. 
The  integration  itself  is  quite  elementary.  However,  the  domain  of  integration  is  quite  complicated 
and  is  bounded  in  general  by  two  straight  lines  and  two  parts  of  a  hyperbola,  if  one  introduces  the 
variable  2  mv^p  for  p.  A  number  of  cases  has  to  be  distinguished  and  the  integration  carried  out 
separately  over  different  parts  of  the  domain  of  integration  before  one  obtains  (5)  and  (5a). 

3.  The  kernel  S  is  greatly  simplified  in  the  case  m  =  1,  i.e.  when  the  moderator  is  hydrogen  gas. 
We  then  have 


yT S(x,xx)  =  4I(x)e  2  (x2  -  xl2) 

2  2 

y7fS(x,x1)  =  4I(xx)e*(xl  x  ) 


(x  <  xx)3 


(x  >  xx), 


ViT V(x)  =  (2x  +  1/x)  I(x)  +  e' 


and  the  integral  equation  becomes 


x 

2x2  J  dxx  I(xx)e  ^xl2  v(xx)  +  4I(x)  eSx  J  e'5Xi2  k(xx)  dxx 


=  (2x  +  )  I(x)  +  e'x2  +  ff  rj  «(x). 

The  S(x,xx)  of  (12)  has  the  form  of  a  Green’s  function  of  an  ordinary  second  order  differential 
equation.  The  integral  operator  of  (14)  is  therefore  the  reciprocal  of  a  second  order  differential 
operator.  If  one  applies  this  operator  to  equation  (14),  it  will  itself  be  transformed  into  a  second 
order  differential  equation. 

Without  using  the  theory  of  the  Green’s  function,  one  can  proceed  as  follows.  Let  us  find  a 
second  order  differential  operator 

h2  H 

L1  =  -  +  a(x)  —  +  b(x) 


(15) 
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which  gives  zero  if  applied  either  to'I(x)  e2X  or  to  e“2  x  .  This  condition  gives  two  linear  equations 
for  the  two  unknowns  a  and  b.  One  obtains 


-  2I(x) 

’x2  +  2xl(x) 


e“x  +  2xl(x) 


In  order  to  bring  into  a  self-adjoint  form,  one  may  multiply  it  with  P~*  where 


P  =  e~x  +  2xl(x). 


This  gives 


d  1  d  e"x 
dx  P  dx  +  p2 


If  one  applies  either  or  L  to  (14),  the  integral  will  give  rise  to  two  types  of  terms.  Differentiation 
of  the  function  before  the  integral  sign  gives  zero.  The  other  terms  will  contain  no  integral.  Hence  the 
integral  equation  reduces  to  a  second  order  differential  equation,  the  second  derivative  coming  from 
applying  Lj  or  L  to  the  right  side  of  (14).  By  virtue  of  these  remarks,  one  sees  that  any  solution  of 
the  integral  equation  (14)  also  satisfies  the  differential  equation 


(V  +  r  )  ^(x)  + 


+  r )  -  4/ 7r 2 


W(x)  =  x2/ P  -  e”x/P2.  (18) 

Since  S(0,X2)  =  0,  we  also  have  y(0)  =  0.  Conversely,  any  solution  y(x)  of  the  differential  equation 
(17)  such  that  v(0)  =  0  is  a  solution  of  the  integral  equation  (14).  The  other  solutions  of  (17)  are,  how¬ 
ever,  not  solutions  of  (14).  The  two  solutions  of  (17)  behave,  at  x  =  0,  like  a  constant  and  x  itself.  The 
former  solution  would  give  an  N(v)  that  goes  to  zero  as  v,  the  latter — the  one  which  we  have  to  use  — 
gives  an  N(v)  that  goes  at  v  =  0  to  zero  as  v^,  which  is  the  same  way  as  the  Maxwell  distribution 
M(v)  goes  to  zero. 

4.  It  would  be  possible  to  discuss  (17)  in  a  general  way  and  to  use  it  directly  to  obtain  numerical 
results.  However,  we  found  it  more  expedient  to  proceed  as  follows.  We  place 


so  that  we  have 


m(x)  =  [v(x)  +  r  j  i/(x), 
(/// P)‘  +^W(x)-4ir‘5  |V(x)  +  /j  ii  =  0, 


M(0)  =  0, 


observing  that  for  x  —  0,  V(x)  >  0,  so  that  V(x)  +  r  ^  0.  We  reduce  the  second  order  linear  differential 
equation  (19a)  to  a  Riccati  equation  by  making  the  substitution 


y  =  fi'/hP, 


(20) 
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which  leads  to 


y  =  W(x)  -  4?T5  [  V  +  r]  -1  -  Py2. 


The  boundary  condition  f±( 0)  =  0,  together  with  |x(x)^  0,  implies  that  we  have 


im  jy(x)  - x-ij. 


Conversely  we  obtain  a  solution  fi(x)  from  a  solution  y(x)  of  (20a)  by  setting 


ju(x)  =  jot'  (0)x  exp  J  [p(t)y(t)  - 1'1  ]  dt. 


The  solution  y(x)  of  (20a)  and  (20b)  may  be  written  in  the  form 
1  4(1  +  6 )x  (103  +  380  8  +  364  6  2)x3 

- -ravor  *  90(1,2i)2  PH 

-  (1163  +  6666  8  +  13200  8  2  +  8864  63)x5 
945(1  +  2  8  )3 

wh  re  6  =  -i —  yW  F  .  We  use  this  series  to  compute  the  values  of  y(x)  for  8  =  0.1,  0.2  and  1  and  for 
x=  0.1,  0.2,  0.3,  0.4,  and  0.5.  With  these  values  of  y  we  compute  successively  /i,  v  and  N,  normalizing 
the  functions  so  that  N"  (0)  =  2. 

With  the  value  of  y(|)  thus  obtained  we  continue  the  function  y(x)  by  numerical  integration  of 
(  :0a).  In  order  to  decide  how  far  to  carry  the  numerical  integration,  we  now  investigate  the  behavior 
of  the  solutions  y  of  (20a)  as  x  approaches  a?  We  use  the  following  expansion  of  I(x) 

I(x)  =  Jf“  --^sr  (1~^i  +-*)* 

2 

Neglecting  terms  involving  e"x  we  see  that 

P(x)  »  yix,  V(x)  »  X  +  ,  W(x)  »  —  (22) 

0T 


Hence  (20a)  becomes 


ym v  =  x  - 


+  4+  i  X 


This  has  two  solutions  which  are  finite  at  ee  and  which  can  be  expanded  in  terms  of  inverse  powers 
of  x: 


8] 
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_  _  ,  2  2  r  (2 r2  +  3)  (2r3  +  5 n 

K  1  2  3  4  5 

xxx  x° 


(23a) 


-!  +  -5 - r- 


2  r  (2r2-i)  (2  r3-T) 


x3  x3 


(23b) 


These  two  series  give  the  asymptotic  expansions  of  two  linearly  independent  solutions  of  (19a): 

/ 1 1  =  exp  j  y1  Pdx  «*exp  (lx2  -  2  In  x  -  2T/x  +  ....)  ,  (24a) 

/x2  =  exp  j  y2  Pdx  **  exp  (-fx2  +  2  In  x  +  2T/x  -  ....)  ,  (24b) 

For  the  asymptotic  behavior  of  the  two  linearly  independent  solutions  of  (17),  this  gives  because  of 
(19) 


1  e(|x2-2  r/x) 

t3  l+Tx-l  +  fx-2 

_(-|x2  +  2  r/x) 


v2*x- 


1  +r  x-i  +|x-2 


and  for  N 


Ni 


1  e 


-2  r/x 

i  +r/x... 


(25a) 


(25b) 


(26a) 


2  e~x2  +  2r/x 
2  iTr7x... 


(26b) 


It  is  evident  from  physical  considerations  that  if  y  >  0  the  slowing  down  at  high  energies,  Nv^,  cannot 
be  zero.  Hence  in  the  asymptotic  expansion  of  N  at  x  =  oo  the  coefficient  of  cannot  be  zero.  Hence 
the  coefficient  of  jjl^  in  the  asymptotic  expansion  of  p,  also  is  not  ero.  As  a  result,  gives  the 
asymptotic  behavior  of  /i  because  for  large  x  we  have  -c  \iy  It  follows  that  it  is  necessary  to 
carry  out  our  numerical  integration  only  to  a  point  Xq  at  which  the  series  (23a)  gives  the  same  value 
as  the  numerical  integration.  For  x  ^xQ  the  series  (23a)  may  be  used  to  obtain  values  of  y.  In  prac¬ 
tice  it  turns  out  that  Xq  =  5  is  a  suitable  choice. 

Having  computed  y  in  this  manner  for  0.5  we  can  compute  /i(x)  by  the  formula 

x 

|i(x)  =  M(i)  exp  j  P(t)  y(t)  dt 

1 

2 


and  then  we  can  compute  i/(x)  and  N(x). 

The  asymptotic  form  of  N  at  large  x  or  large  v  is  given  by  (26a).  Replacing  e~2^Vx  in  this  by 
(1  +  r  /x)~2 — which  is  within  the  accuracy  of  the  expression  in  the  denominator — we  obtain,  rein¬ 
troducing  4  s  and  the  life  time  r  of  the  neutrons, 
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N(v)  dv 


(Nffc  T  v  +  1)' 


where  C  is  a  constant.  This  shows  that  the  deviation  from  the  dv/v2  law  arises  at  high  energies  not 
from  the  positive  temperature  of  the  atoms  but  from  the  absorption.  It  causes  N(v)  to  become  lower 
than  the  dv/v2  law  would  indicate.  At  lower  energies  the  Figures  show  that  N(v)  becomes  higher  than 
1/v2.  This  is  the  result  of  the  finite  temperature  of  the  atoms.  For  T  =  0,  (27)  is  correct  for  all  v. 
The  constant  C  occurring  in  (27)  is  of  some  importance  and  may  be  evaluated  directly  from  the  values 
of  N(x)  previously  calculated.  If  this  be  done  it  is  possible  to  obtain  a  check  on  the  computations  used 
to  find  N(x)  by  calculating  C  in  another  manner.  The  number  of  neutrons  which  are  slowed  down  to 
velocities  below  v0  is  given,  when  vq  is  very  large,  by  the  formula 


r  vN(v)vog_dv^  lim  v2n(v) 

S  J  2  « 


2N2  a  2t3 
s 


Since  vQ  is  large,  this  is  approximately  equal  to  the  number  of  neutrons  that  are  absorbed  at  velocities 
below  v0  and  in  the  limit  as  vQ, approaches  infinity,  we  have  equality.  Thus 


N asy  J  N(v)dv 


2N2  a  2  t3 


Hence  the  constant  C  of  the  asymptotic  expression  for  N(v)  can  also  be  obtained  by  integrating  N(v). 
The  integral  was  evaluated  as  follows.  Replace  vbyx  =  0v,  so  that 


N(x)dx 


2N3as3  t3 


The  integral  ]  N(x)  dx  =  A  may  be  computed  by  numerical  integration  of  the  values  of  N(x)  previously 
0  oo 

computed.  To  compute  B  =  j  N(x)  dx,  we  replace  N(x)  by  an  asymptotic  expansion  of  the  form 


N(x)  «  — ^ -  —  + 

2N3  or  3  T3  x2 


(which  may  be  obtained  by  making  more  precise  the  reasoning  which  leads  to  26a)  and  integrating 
directly.  In  this  manner  we  get  a  linear  equation 


A  +  B  = 


2N3cts3  t3 


which  may  be  solved  for  C.  The  constant  C  was  computed  in  this  manner  for  the  three  values  of  y 
previously  mentioned,  and  it  was  found  that  they  agreed  satisfactorily  (about  1%)  with  the  values  of  C 
computed  directly. 
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Our  numerical  results  are  summarized  in  the  attached  Figure  which  shows  the  graphs  of 

9  2  1 

M(x)  =  x*e“x  (y  =  0)  and  N(x)  for  5  =  —  y7  0  y  -  0.1,  0.2  and  1  respectively.  This  corresponds  to 

ffa/ *s  valnes  0.226,  0.452,  2.257  at  kT  energy.  In  the  three  latter  cases,  the  straight  line  graph 

of  ~1T  x2n(x)  is  aIso  included,  so  that  the  deviation  from  the  dv/v2  law  may  be  visualized. 

2  x~ 00 


Figure  1.  Graphs  showing  N(x)  for  0y  =  0,  0.226,  0.451,  and  2.257,  and  its  asymptotic  behavior  for 
large  X  for  fiy  -  0.226,  0.451,  and  2.257.  x  is  proportional  to  the  velocity,  X  =  1  corresponding  to  kT. 
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